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ABSTRACT 


The radar cross section (RCS) of various geometries with 
surface resistivity is synthesized using a physical optics 
approximation. Using the RCS, a set of equations are 
developed for reconstructing the resistivity across the 
surface. The reconstructed resistivity is then compared with 
the original resistivity and the RCS recalculated using the 
synthesized resistivity to validate the accuracy of the 
synthesis procedure. 
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I. INTRODUCTION 


During the middle 1950's, the idea of making platforms 
unobservable or less observable to radar was initiated by Bill 
Bahret of Wright Air Development Center. [Ref. 1] The initial 
motivation for low-observability came from the need to 
increase screening provided by a jammer. [Ref. 1] 
Subsequently, a great deal of research in the area of radar 
cross section (RCS) began. 

There are four common methods for reducing the RCS of a 
given platform. They are [Ref. 2]: 

1. shaping, 

2. surface material selection, 

3. active cancellation, and 

4. passive cancellation. 

Shaping to reduce RCS requires reducing the surface area that 
is perpendicular to the radar source. Surface material 
selection requires the addition of lossy materials to the 
surface of the target. Active cancellation is a method of 
reproducing the signal 180° out of phase from the radar. 
Passive cancellation is the addition of a secondary scatterer 
to cancel the reflected field from the primary target. 

In this thesis, a synthesis method that can be used for 
surface material selection is addressed. That is, given the 
monostatic RCS of a platform, a set of equations are presented 
which can be used to choose the required material. Selecting 
a surface material that reduces the RCS of a platform requires 
the use of a radar absorbing material (RAM). Other types of 





materials used in the reduction of RCS include resistive and 
conductive films, described in terms of ohms per square (Q/D) 
or mhos per square (U/D), respectively. 

Metal surfaces such as those used for aircraft and ships 
closely approximate a perfect electric conductor (PEC). The 
addition of RAM to a surface causes it to have an inductive or 
capacitive surface impedance rather than that of a PEC. 
Impedance surfaces can absorb radar frequency (RF) energy 
thereby causing a reduction in the scattered energy. 

An impedance surface supports both equivalent magnetic 
and electric currents. However, there are two special cases 
of impedance surfaces that only support either equivalent 
magnetic or electric currents: a resistive surface does not 
support equivalent magnetic currents, and a conductive sheet 
does not support electric currents. The resistive surface is 
the focus of this thesis, because its use greatly simplifies 
the equations used in the synthesis procedure. 

A synthesis procedure can be developed based on the 
desired RCS for a particular target. The surface is separated 
into triangular subsections, or patches. The surface 
impedance across each subsection is assumed to be constant but 
can vary from patch to patch. A matrix equation is derived 
consisting of a vector of the known scattered field, a target 
surface scattering matrix, and a surface impedance vector. By 
solving the matrix equation, the impedance vector is obtained 
and the material properties of the surface can be 
determined. 

In Chapter II, the theoretical background required to 
evaluate the RCS of a surface is presented. Surface 
impedances, reflection and transmission coefficients and the 
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physical optics approximation are examined. 

In Chapter III, the synthesis procedure is developed. 
The procedure is applied to doubly curved surfaces such as 
spheres or portions of spheres. 

In Chapter IV, the results of the synthesis procedure are 
compared to the exact solution. Finally, in Chapter V, an 
analysis of the results is given as well as recommendations 
for future work in this area. 







n. THEORETICAL BACKGROUND 


In th-is chapter, the basic equations of electromagnetic 
scattering are presented. 

A. RCS FROM THE SCATTERED FIELD 

The RCS of a target can be determined once its scattered 
field is known [Ref. 3] 

lim (2.1) 

The scattered field arises from the currents induced by the 
plane wave on the target surface. The scattered electric 
field far from the target is given in terms of the electric 
(J) and equivalent magnetic (M) currents by [Ref. 4] 

. ( 2 . 2 ) 

ke 

E*=j-!^(Le-T,N^) . 


Ng= Jj (J^cos6cos(|)+JyCos6sin(J)-J^sin0)e-’^' “‘^ds 
N ^=f J(-J^sm<l)+JyCos<|>)ej‘'''“*ds' , 


(2.3) 

(2.4) 

(2.5) 











B. 


SURFACE IMPEDANCES 


Traditionally, most aircraft and ships were made of 
materials which are considered to be perfect conductors. 
Conducting targets generally have a high RCS relative to those 
of other materials. In RCS reduction (RCSR), the selection of 
materials for a target can reduce the scattered field. One 
method is to coat a perfectly conducting target with RAM. 
This yields a heavy target and is not an effective way to 
reduce RCS. Another approach is to use composite materials 
with some loss introduced. The latter approach is more 
efficient, since composites are frequently used in the 
manufacture of aircraft and ships. Composites can be 
characterized by their complex index of refraction, n,, which 
is related to the permittivity and permeability of the 
material 


( 2 . 8 ) 


where e, is the permittivity of the material and n, is the 
permeability of the material. The permittivity and 
permeability are usually complex numbers and it is the 
imaginary parts which contribute to losses. [Ref. 2] Thin 
layers of lossy material can be applied to a dielectric to 
cause the surface to be resistive and, no magnetic currents 
will flow on the surface. 

If the target consists of composite materials then there 
is some penetration of the field into the target. The 
integrals in equations (2.4) and (2.5) become volume integrals 
because the currents exist below the surface. This is very 







inconvenient computationally, and the problem can be avoided 
using the surface impedance approximation, also known as the 
Leontovich boundary condition 





(2.9) 


E^„ and are the total tangential fields on the outer 
surface of the target and t), is its surface impedance in ohms 
(ri,= 0 if the surface is a PEC) . [Ref. 3] Thus the original 
target surface can be replaced by a surface with impedance v), 
and the same scattered fields result. Note that r|, depends on 
the location on the surface and the incidence direction. 

Using equation (2.9), a set of vector equations can be 
derived 


E-a (fit • 


(2, 

.10) 


_ M 

(2, 

.11) 

H-A (fi ■ 

. H)=—i . 


For a resistive surface, this simplifies to 

fi X (£i X E)=-r7, , (2.12) 

fi X (£i X H)=0 , (2.13) 


where R, is the resistivity of the material. 



C. REFLECTION AND TRANSMISSION COEFFICIENTS 

The reflection coefficients quantify the amount of 
reflection and transmission of an electromagnetic wave that 
occurs at an interface between two different impedances. For 
a PEC, the incident wave is totally reflected, yielding a 
reflection coefficient of -1. The reflection (P) , 

absorption (A) , and transmission (T) coefficients must obey the 
conservation of energy as shown below [Ref. 2]: 

|r|^+|A|^-.-lTl^=l . (2.14) 

For a lossless material, A=0. 

The transmission and reflection coefficients are easily 
derived with the assumption that the surface is locally flat. 
Two principle polarizations must be considered: 1) transverse 
magnetic to the z axis (TM) and 2) transverse electric to the 
z axis (TE). An arbitrary polarized wave can be handled as a 
combination of the two principal polarizations. 

1. TM Polarization 

For TM polarization, the incident wave is 

Ei=0Eee'''‘’' (2.15) 

where kj is the incident wave vector. For a TM polarized 
incident wave, the reflection coefficient for a resistive film 


[Ref. 






-ti|jCos0j 

2Rj+t1j^cos0j 


(2.16) 


where 0j is the angle of incidence with respect to the z axis. 
If no absorption occurs, as is the case with an infinitely 
thin sheet, the transmission coefficient is 

T (2.17) 

™ 2R^+ii^cos0; 

The scattered field then becomes 

E.=0r.j.^E,e'-’‘‘' (2.18) 

where k is the wave vector in the direction of observation. 

2. TE Polarization 

For TE polarization, the incident wave is 

ij=(|>E^e’^‘‘‘' . (2.19) 

For a TE polarized incident wave, the reflection coefficient 
is [Ref. 5] 


( 2 . 20 ) 


With no absorption, the transmission coefficient is given by 
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(2.21) 


2R^cos0j 

2R^cos0j+ti^ 


and the scattered field becomes 


( 2 . 22 ) 


D. PHYSICAL OPTICS APPROXIMATION 

The field incident on an illuminated object causes 
currents to flow on its surface. To solve exactly for the 
scattered electric field, these currents must be known. It is 
a lengthy and complicated process to solve for the currents. 
However, the incident field and the boundary conditions are 
known. Therefore an approximation can be made using the known 
incident field. The standard physical optics approximation is 
given by [Ref. 3]: 


J,= 2 fi X Hj . (2.23) 

Equation (2.23) holds for the portion of the object which 
is illuminated. Physical optics approximation assumes the 
current along the surface drops to zero on areas of the object 
which are shaded from the illuminated field. In actuality, 
currents do not drop abruptly to zero. The currents slowly 
decay towards zero, but for electrically large bodies the 
error is negligible. 
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The physical optics approximation does not account for 
waves traveling along the surface of the object. Therefore, 
there is a decrease in RCS as the incident angle moves towards 
surface grazing and the sidelobes of the RCS pattern are lower 
than the actual values. The specular lobe level of the RCS 
pattern for flat or gently curved surfaces is generally within 
a few percent of the actual value. 



m. FORMULATION AND SOLUTION 


The analysis of scattering from an arbitrary finite 
doubly curved surface requires the use of triangular 
subsections. The physical optics approximation and boundary 
conditions are used to obtain the synthesis equations. The 
synthesis equations are then verified by comparing with 
calculated fields. The geometries considered for this 
calculation are a flat disk and a spherical cap. The 
geometries are shown in Figures 2 and 3. 

The incident wave is assumed to be a planar wave with 
arbitrary polarization. Both TM and TE components exist. The 
incident field is given by 

V(0eMe;) (3-1) 


-jkj«r=jk(xsm0;Cos<J>i+ysm0jsin<{)j+zcos0;) . (3.2) 

Two cases are considered for the formulation and 
solution: bistatic and monostatic. Bistatic is the most 
general case. The transmit and receive antenna for the radar 
are not located at the same position. Monostatic is the case 
in which the transmit and receive antenna are co-located. 
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For each case, a flat disk is considered first. The disk 
is a flat surface with a curved perimeter that requires the 
use of triangular subsections to model. Next, the more 
complex shape of a doubly curved surface is solved. The 
solution for a doubly curved surface is more complex because 
a transformation of coordinates is required for each 
triangular subsection. This transformation gives the surface 
normal in the transformed z axis direction. Also, for the 
doubly curved surface, the incident wave polarization in the 
individual triangle coordinate systems is usually not pure TM 
or TE, even if the polarization is purely TM or TE in the 
global coordinate system. 

A. BISTATIC CASE 
1. Flat Disk 

For a flat disk in the x-y plane, the normal is always in 
the z axis direction. Therefore, the boundary condition from 
equation (2.9) can be written as: 

i X (i X E)=-R,J. (3.3) 

where the electric field in the above equation is given by 

E,= (0T^^E^$T„E;)e ^ . 4) 


The triple products in equation (3.3) becomes 
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3; X (4 X E)=[(-ftcos0jSm({).+^cos0jCOS<t>i)T^^Ea'+ 

(-icos0jSin<{)j-^cos<{)jCOs0j)T^E^‘]6 . (3.5) 

By substituting and T.[j from equations (2.25) and (2.29) 
into (3.5) and dividing by -R,, equation (3.3) can be solved 
for the current along the surface 


2Rj+ili^cos0j 

2(-tsm<t)jCOS0j-^cos0jSiint)j) 


(3.6) 


The two principal polarizations are considered separately. 


a. TM Polarization 

For TM polarization, the (|) component of the incident 
field is zero. Therefore, J, in equation (3.6) simplifies to: 


- 2{-tcos0;Cos(t)j+^cos0;Sin4);) 
‘ 2R,+ti^cos0j 


For use in equation (2.4), equation (3.7) is separated into x 
and y components 


-2cos0jCos<{>j i .j£_., (3.8) 

^ 2R,+tl^cos0j ® 

2cos0jSin(J); j .jj., 

2R,+Tij^cos0j ® (3.9) 
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Equation (2.4) simplifies to; 


rr 2COS0.COS(().COS0COS<{> j 

Ne = i J [-^^^ 

2R.+ti^cost0i 

2cos0.sin<t).co50sin(t> j , .jf,., ^ ^ 

-‘ dx'dy 

2R,+’l„cos0j 


Portions of equation (3.10) are not dependant on surface 
location and can be removed from the integral. Equation 
(3.10) then reduces to 


N,=2EgCos0jCos0cos(<()j+<t)) 

,-jh-^jkr'<=»•♦ , , 

X I /-dx dy . 

•^•'•2R.+ti„cos0i 


The portion of equation (3.11) inside of the integral can be 
denoted as I for simplicity 


■=//. 


-dx'dy' . 

'2R,+ti^cos0; 


(3.12) 


The resistivity of the material is restricted to be 
constant on each triangle, but can be a can be a function of 
location on the surface of the object. Therefore, it must 
remain within the integral. Thus the resistivity function can 
be defined as a series of pulses with unknown coefficients 
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(3.13) 


R.(’^c.yc)=E 

where c is the triangle number and N is the total number of 
triangles on the surface of the object. The coefficients 
are the unlcnown expansion coefficients and f^. are the step 
functions 

(l. 

t=\ (3.14) 

Ip, otherwise 

where is the area of triangle c. Substituting equation 
(3.13) into (3.12) gives 

1 = j -dx'dy . 

(3.15) 

Brealcing the integral into a sum of integrals over the 
individual triangles gives 




—dx'dy's^I^ 
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The c'*’ term of the sum is 


2r^ + r|^cos0i 


(3.17) 


The denominator of equation (3.17) can be defined as 

1 (3.18) 

° 2r_^+Ti^cos0j 

The integrals over the individual triangles can be reduced to 
closed form as described by Moreira and Prata. [Ref. 6] In 
general, 

e'’^’C C (3.19) 

I^=2S„e^^-[-^-2-5-] . 

Dj,(D^-Dp D^(D^-Dj,) 

In equation (3.19), the C and D coefficients are given by the 
following equations: 




(3.20) 

(3.21) 

(3.22) 

(3.23) 
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Case 1, lDpl<Lt and |D,|iLt: 


+e^'CpG(n,-D^)] . 


(3.24) 


Case 2, | Dp |<Lt and |D,|<Lt; 


•EE^f^Cp. 

n.o m-o (n+m + 1)! 


(3.25) 


Case 3, |Dp|>Lt and lD,|<Lt: 




—G(n + 1,-Dp)] . 
n + 1 


(3.26) 


Case 4, IDpULt, ID^IiLt, and |Dp-D,|<Lt; 


UCDp-Dp)]" e^’Cp 

Io=2Sp—E- - -3—[-CpG(n.Dp) +-2] . 

“ ‘jD n! ° ^ (n + l) 


(3.27) 


In the above equations, n and m are the series terms, Lt is 
a small number that determines the region of applicability of 
the series approximation, and G is a function of the form: 


G(n,w)=/‘s “ej"'ds . 


(3.28) 
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For the calculations in this thesis, values of Lt=.005 and 
m=n=2 were sufficient. 

Substitution of equations (3.16) and (3.19) into (3.10) 


Ng=2EaCos0jCos0cos((t>j-4>)53 


(2.2), the expression for the co-polarized 


—EgCos0jCos0cos(4>;-<{>)5^ 


Finally, using I in place of the 


To calculate the cross-polarized scattered field, 
equation (2.5) is used and the same procedure employed to 
obtain the co-polarized RCS is followed. Equation (2.5) 


Equation (2.3) for the scattered field become 
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and the cross-polarized RCS is 




lcos0jSin(<t)-<t>;)I|^ 


b. TE Polarization 

For TE polarization, the 0 component of the incident 
field is zero. Therefore, J, in equation (3.6) simplifies to 

_ 2(-tsm<t)jCOs0;-^cos0.sin({)j) ; ^ 35 ) 

‘ 2R^cos0j + Ti^ * 

For use in equation (2.4), (3.35) is separated into x and y 

components. 

^ _ -2(sm»,cos0,) ^,^,jg,., ^ (3.36) 

^ 2R^cos0j+no * 

^ -2(cos6,gm({>,) ^ (3.37) 

2R^cos0j+Ti^ * 

Equation (2.4) simplifies to 

-- jii'•<>•♦ -j»r' (3 381 

N-=-2(cos0cos0.sin(<t)-4).) I ( -dx'dy' . 

•'•'• 2 Rcos 0 j+n„ 
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The equation for I corresponding to equation (3.14) becomes 


r r - jir'oo'tp 

1=1 ^^-djcMy' . 

•' •' •2RjCOS0;+ti^ 


Equation (3.19) for I,, remains the same. However, the 
denominator is not the same and therefore, a new variable b. 


The equation for the cross-polarized scattered field becomes 


the cross-polarized RCS : 


The scattered field for the co-polarized directii 
derived in a similar manner yielding 
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c. Arbitrary Polarization 

The most general incident polarization is a 
combination of TM and TE components. The total scattered 
field in the 0 direction is 


Ee 




(3.45) 


—Ejlcos0j[cos0cos(<t)-({)j)+sin(4)-<{)j)] . 


(3.46) 


And the total scattered field in the (j) direction is 


(3.47) 


or. 


—E^‘lcos0j[cos0sin(<J)-(|)j)+cos(<()-<()j)] . 


(3.48) 




Equations (3.47) and (3.48) 


be put into a matrix form. 



(3.49) 


where 


—EjCos0jCos0cos(<J)-<J)j) , 


(3.50) 


—E^cos0jCOS0sin(<{)-<{)j) , 


(3.51) 


A 


2jir 


Eg‘cos0;sin(<{)-4)j) , 


(3.52) 


-jkT1„e- 


-E^'cos0jCos((t>-(J)j) , 


(3.53) 




(3.54) 
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and 



The resistivity of each triangle is embedded in and via 
equations (3.54) and (3.55). 


2. Doubly Curved Surfaces 

For a doubly curved surface, the surface normal vector of 
each triangular subsection can be different. The equations 
(3.19) through (3.27) still hold if the local triangle 
coordinates are used. Therefore a transformation of the 
coordinates must be performed for each triangular facet on the 
surface from the global system (x,y,z) to the local system 
(x",y",z") as shown in Figure 5. The z" axis is in the same 
direction as the facet normal n^. After the transformation of 
coordinates, the scattering calculations are made for each 
triangle, the fields summed over all triangles, and then a 
transformation is made bacic to the original (global) 
coordinates. 

The incident field is first converted to global Cartesian 
coordinates [Ref. 4] 

cos(t>sinO co50cos<t> -Ein(t> 
sm6sin4> cosBsinct) cob4> 
cos6 -sinO 0 
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Figure 5. Transformation of coordinates from a global (x,y,z) to a local (x",y",z") 
system for triangle c. 

lext, a rotation about the z axis is performed to give a 
;oordinate system l.x\y‘,z') . Then a rotation about the y' axis 
lone to give the local- coordinate system of triangle c (x",y' 







The angle is the angle of rotation about the z axis, and the 
angle is the angle of rotation about the y' axis, where 



and 


p„=cos-*(ft^ X 4) . (3.59) 

The incident field is also transformed to the local 
coordinate system using the same matrices as in equation 
(3.57). After all transformations are complete, equations 
(3.45) through (3.55) are applied. Once the scattered field 
for each triangle is calculated, it is transformed bade to the 
global coordinate system and the fields from all triangles are 
summed vectorially. Equations (3.31), (3.34), (3.42), and 

(3.45) are then used to calculate the total RCS of the 
surface. 


B. MONOSTATIC CASE 

For a monostatic radar, <!);=<}) and 0i=0 and the scattered 
fields reduce to 


E,'= 


2iir 


Eglcos^e 


(3.60) 
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. -jkti e'^ : (3.61) 

E>-5-E.Icose . 

* 2itr * 

In equations (3.60) and (3.61), I is still given by equations 
(3.15) and (3.38). Recall that and are dependant on 0j 
and location on the surface. In order to synthesize an 
equation for resistivity, and cannot be dependant on 0;. 
This implies that the surface composition must change for each 
incidence angle (ie., an adaptive target). One approximation 
that eliminates the dependence on 0; is to consider cos0j=l. 
In this case. 


and equations (3.60) and (3.61) reduce to 


(3.63) 


E, 


2itr 


-Ski 

—e;i . 


C. SYNTHESIS EQUATIONS 

In applying the synthesis procedure, only the monostatic 
case is considered. The scattered field at various angles is 
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To generate field values for the left sides of equations 
(3.63) and (3.64), DCS is divided into squares of area du by 
dv. There are Nu divisions along -the u axis and Nv divisions 
along the v axis. The center of each square is the sample 
point for calculating the scattered field. This gives a total 
of Nu by Nv samples of the scattered field. However, only the 
points within the visible region are valid for use in the 
synthesis equations. The nximber of points inside the visible 
region is defined to be N. The angles 0 and (f) for use in 
calculating the scattered field are found from 

e = sin'‘(u^+v^) (3.65) 


The unlcnowns in the synthesis problem are the expansion 
coefficients, a^. There is one value of for each triangular 
subsection. Let a be a vector of N elements representing N 
un)cnowns. In order for the equation to be solved uniquely, N 
un)cnowns require N equations and, therefore, the number of 
samples of the scattered field must also be N. A matrix 
equation follows from equations (3.15), (3.60), and (3.61); 
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En=E ■(e„,4>^) = E '(u^.v^) , 


and c,n=l,2, . . .,N. 

The X matrix is square and can be inverted and used to 
solve for the a matrix 


X„ X,j 

x„ 


Then using equation (3.62), each element in the a matrix is 
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used to solve for the resistivity over each triangle. 


Equation (3.72) represents the synthesis equation for 
calculating the resistivity of an object given the required 
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rV. COMPUTER IMPLEMENTATION AND DATA ANALYSIS 


This chapter discusses the computer coding of the 
equations derived from Chapters II and III. The first 
sections cover the calculation of the RCS, while the second 
half addresses the synthesis calculation. 

A. RCS FROM THE SCATTERED FIELD 

Using equation (2.1), the RCS can be calculated from a 
scattered field. The incident field is assumed TM polarized 
and the monostatic RCS is calculated for the scattered field 
obtained using equations (3.60) and (3.61). A listing of the 
program POTEST is included in the appendix. It should be 
noted that POTEST calls several subroutines which are also 
included in the appendix. 

A lA X IX flat PEC plate was used as the initial test 
case. The patch model of the plate is shown in Figure 7. In 
Figure 8, the approximate (PO) results are compared to the 
method of moments (exact) results. The specular lobe was 
within a few dB, but the sidelobes are slightly larger in the 
exact results. This is due to the fact that the approximation 
used in equation (3.3) neglects any traveling waves. 

The second test was to place some resistivity on the 
surface of the plate. First, a constant resistivity of 377 
ohms was used. A comparison of the PO and MOM results is seen 
in Figure 9. Note that with a resistance on the surface, the 
sidelobe differences decrease. 
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The final test was to increase the resistivity linearly 
from 0 to 377 Q/D from the center of the plate to the edges. 
This was done using the equation 

R.=377(x,^+yc^) {4.1) 

where andy^ are the midpoints of triangle c. Equation (4.1) 
is for a flat surface laying in the x-y plane, centered at the 
origin. The results are shown in Figure 10. 

The next step was to look at a flat disk and a spherical 
cap with and without resistivity to establish the accuracy of 
PO relative to MOM. The models used for a flat disk and a 
spherical cap are shown in Figures 2 and 3, respectively. The 
results for a flat disk are shown in Figures 11 through 13. 
Figure 11 is the RCS of a conducting disk. Figure 12 for a 
disk with constant resistivity of 377 Q/D, and Figure 13 is 
for a disk with resistivity given by equation (4.1). Figures 
14 through 16 are the results for a spherical cap. Figure 14 
shows the RCS for a PEC cap. Figure 15 for constant 
resistivity of 377 Q/D across the cap, and Figure 16 for 
linear resistivity using 

R.=377^ (4.2) 

where is the arclength from the pole to the center of 
triangle c and, is the total arclength from the pole to the 


base 


the cap. 




Theta in Degrees 


Figure 10. RCS of a flat plate with surface resistivity given by equation (4.1)(TM polarized 
incident wave, monostatic). 



































B. 


RCS SYNTHESIS 

The synthesis computer code is similar to the program 
POTEST with changes to utilize DCS samples rather than (]) and 
0. The increments du and dv were calculated so that N would 
equal the number of triangles on the surface of the object. 
A short program, NBOX, was used to insure this. After the 
number of boxes and triangles were matched, the program P0TC02 
was used to calculate the field scattered from the object at 
the required u and v points. The scattered field was then 
used in SYNTH to calculate the resistance on each triangle 
based on the given field. Listings of all of these programs 
are included in the appendix. 

First, a flat conducting plate was again used as a test 
case. For this calculation, N=80, du=dv=0.2. The synthesized 
resistivity was expected to be 0 0/D across all triangles, and 
this was indeed the result. Next, the plate was given a 
constant resistivity of 377 0/D. The synthesized resistivity 
exactly matched the expected original values. Finally, the 
plate was given a linear resistivity across the triangles. 
The synthesized resistivity was within less than one percent 
of the expected result. A comparison plot of the original 
resistivity and the synthesized is shown in Figure 17 on the 
following page. Note that adjacent triangles have different 
resistivity curve oscillates. 
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The flat disk was next used to calculate the synthesized 
resistivity. For this calculation, N=156, du=dv=0.1429. With 
a constant resistivity of 0 fi/D and 377 Q/D, the synthesized 
resistivity was within one percent of the original value. 
This is shown in Figures 18 and 19 respectively. With a 
linear resistivity given by equation (4.1), the synthesized 
resistivity was exactly the same as the original resistivity, 
as shown in Figure 20. 

Finally, resistivity synthesis was conducted for a 
spherical cap. For this calculation, N=156, du=dv=0.1429. 
The synthesized resistivities for cases 0 Q/D and 377 fi/D 
were within one percent of the original resistivities. This 
is shown in Figures 21 and 22. The results for a linear 
resistivity given by equation (4.2) were exactly the same as 
the original values. This is shown in Figure 23. 

In order to verify the synthesized resistivity would 
yield the correct RCS, a comparison was made between the 
original RCS and the synthesized RCS. For all three 
geometries, 0 fl/D, 377 D/D, and linear resistivity, the 
synthesized resistivities were used in the MOM code to compute 
the RCS patterns. The results are shown in Figures 24 through 
32. The curves labeled Orig are the PO RCS values used to 
initiate the synthesis. Those labeled Synth were obtained 
using the synthesized resistivity in the MOM code, PATCH. 
These figures are essentially duplicates of Figures 8 through 
16 because the synthesized resistivity was almost identical to 
the starting values. 
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V. CONCLUSIONS 


A procedure for synthesizing the RCS of a target of known 
shape has been presented. The equations were derived for both 
flat and doubly curved surfaces and demonstrated for 
monostatic, and purely resistive surfaces. The equations were 
derived using the physical optics approximation and, 
therefore, traveling waves were neglected. Furthermore, 
multiple reflections were neglected. 

The RCS of each type of surface was first calculated 
using PO and then compared to MOM. Next, the RCS of each 
surface was calculated in DCS so the scattered field could be 
sampled and used in the synthesis equations. The synthesis 
equations were then used to find the resistivity across the 
surface. Finally, the resistivity obtained from the synthesis 
procedure was used to recalculate the RCS and, therefore, 
demonstrated that the synthesis procedure was accurate. The 
synthesized resistivities were in excellent agreement with the 
original values in spite of all of the approximations. This 
is because the chosen resistivity distributions yielded low 
RCS sidelobes and did not support traveling waves. Similarly, 
the resistive surface transmits more of the incident field 
than -it reflects. Therefore, multiple reflections are not 
important and can usually be neglected in the synthesis of low 
RCS patterns. 

A more rigorous solution could be obtained by using the 
MOM as the basis of the synthesis equations. [Ref. 3] This 
would require solving exactly for the current along the 
surface as well as for the resistivity. Future work in this 
area should include a synthesis procedure based on MOM. 
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APPENDIX. COMPUTER CODES 


This appendix contains all of the computer codes used for 
calculations in this thesis. The code POTEST.M calculates 
the RCS of a surface with specified resistivity. The code 
P0TC02.M calculates the RCS of a surface in DCS. The code 
BUILD.M generates the geometry of the surface. The code 
FACT.M is a subroutine used to calculate the factorial of a 
number and is used with POTEST.M and P0TC02.M. The code GF.M 
is a subroutine used to calculate equation (3.28) and is used 
with POTEST.M and P0TC02.M. The code RES.M is a subroutine 
used in the calculation of linear resistivity across a surface 
using equations (4.1) and (4.2) and is used with POTEST.M and 
P0TC02.M. The code SYNTH.M reconstructs the resistivity using 
the scattered field from P0TC02.M. The code NBOX.M calculates 
the number of partitions of DCS that are within the visible 
region for a given du(=dv) . This was used to verify that the 
number of Icnown scattered field points would be equal to the 
number of triangular subsections on the surface. The code 
POS.M compares the original RCS with the synthesized RCS. 
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% POTEST.M 

% PO scattering from triangles 
clear all 

% illumination flag: =0 external face only 
iflag=l; 

Lt=.001; 

Nt=2 ; 
wave=l; 
bk=2*pi/wave; 
rad=pi/180; 

% Incindent wave polarization 
Et=l+j*0; %TM-z 

Ep=0+j*0; %TE-z 

phi=0; 

phr=phi*rad; 

cp=cos(phr);sp=sin{phr) ; 

% Wave ampltidue at all vertices 

Co=l.; 

build 

ntria=nfaces; 

nvert=nverts; 

x=xpts; 

y=ypts; 

z=zpts; 

title('triangular model of the scattering surface') 

xlabel('X') 

ylabel{'y') 

zlabel('z') 

for i=l:nvert 

text(x(i)-max(x)/20,y(i)-max(x)/20,z(i), num2str(i)) 
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end 

hold off 

dispChit return to continue calculation') 
pause 

% Define position vectors to vertices 
for i=l:nvert; 

r(i,:)=[x(i) y(i) z(i)]; 

end 

%Use res.m program to calculate resistivity 
res 

% Get resistivity, edge vectors and normals from edge cross 
products 

for i=l:ntria 
rs(i)=0; 

% Resistivity function for a flat surface 
%rs (i) =rss*sqrt (xm (i) ''2+ym (i) ^'2) /377 ; 

% Resistivity function for a doubly curved surface 
%rs (i) =rss*arc (i) 12 , 11 ; 

A=r (vind {i, 2) , :) -r (vind (i, 1) , ; ) ; 

B=r(vind(i,3),:)-r(vind(i,1),:); 

C=r(vind(i,3),:)-r(vind (i,2),:); 

N(i, :)=cross(A,B); 

% check for proper direction of normal vector 
if N(i,3)<0 

N(i,:)=-N(i,:); 

end 

% Edge lengths for triangle "i" 
d (i,1)=norm(A); 
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d(i,2)=norm(B); 
d(i, 3) =norm(C); 
ss=. 5*siHn(d(i, :) ) ; 

Area(i)=sqrt(ss*(ss-d(i,l))*(ss-d(i,2))*(ss-d(i,3) ) ) 
Nn=norm(N(i,:)); 

N{i, :)=N(i, :)/Nn; 

%transforni triangles into local coordinates 
xl=[l 0 0]; 
yl=[0 1 0]; 
zl=[0 0 1]; 

end % end of triangle loop 

% Pattern Loop 
start=0; 
stop=180; 
del=l; 

it=floor((stop-start)/del)+1; 
for i=l:it 

theta(i)=start+(i-1)*del; 
thr=theta (i)*rad; 
st=sin(thr); ct=cos(thr); 
u=st*cp; v=st*sp; w=ct; 
up=ct*cp; vp=ct*sp; wp=-st; 

% Spherical coordinate system unit vectors 
R(i,:)=[u V w3; 

% Change Efield into global rectangular coordinates 
strtra(1,:)=[st*cp ct*cp -sp]; 
strtra(2, :) = [st*sp ct*sp cp]; 
strtra(3, ;) = [ct -st 0] ; 
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sei=[0 Et Ep]; 
ei=strtra*sei' ; 

% correct orientation for normal after 90deg 
if theta(i)==91 
N=-N; 

end 

% Begin loop over triangles 


for m=l:ntria 
t (m, :)=cross(zl,N(m, :) ) ; 
beta(m)=acos(N(m,3)); 

alpha(m)=atan2(t(m,2),t(m,1)+le-6)-pi/2; 

% Test to see if front face is illuminated 
ndotk=N(m,:)*R(i, :)'; 
if ndotk >= 0 & iflag==0 | iflag~=0 
%translation into local coordinates 
%first rotate about z axis 

tral(1,:)=[cos(alpha(m)) sin(alpha(m)) 0]; 

tral (2, :) = [-sin(alpha (m) ) cos(alpha(m)) 0]; 

tral(3,:)=[0 0 1]; 

% Rotation about yp axis 

tra2(1,:)=[cos(beta(m)) 0 -sin(beta(m))]; 

tra2(2,:)=[0 1 0]; 

tra2(3,:)=[ sin(beta(m)) 0 cos(beta(m))]; 
Dl=tral*R(i,:)'; 

D2=tra2*Dl; 
upp=D2(1); 
vpp=D2(2); 
wpp=D2(3); 
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% Find spherical angles in local coordinates 
stpp=sqrt (upp''2+vpp''2) *sign (wpp) ; 
ctpp=sqrt {l-stpp''2) ; 
thpp=acos(ctpp); 
phipp=atan2(vpp,upp+le-10); 

% Phase at the vertices of triangle m; monostatic RCS needs 
"2" 

Dp=2*bk;* ( (x (vind (m, 1) ) -x (vind (m, 3) ) ) *u+ . . . 

(y (vind (m, 1) ) -y (vind (m, 3) ) ) *v+ . . . 

(z (vind (m, 1) ) - z (vind (m, 3) ) ) *w) ; 

Dq=2*bk*((x(vind(m,2))-x(vind(m,3)))*u+... 

(y(vind(m,2))-y(vind(m,3)))*v+... 

(z(vind(m,2))-z(vind(m,3)))*w) ; 

Do=2.*bk* (x (vind (m, 3) ) *u+y (vind (m, 3) ) *v+z (vind (m, 3) ) *w) ; 
%incident field in local rectangular coordinates 
eip=tral*ei; 
eipp=tra2*eip; 

% translation of E-field into local spherical coordinates 
strtra2=[stpp*cos(phipp) ctpp* cos(phipp) -sin(phipp); 

stpp*sin(phipp) ctpp*sin(phipp) cos(phipp); ctpp -stpp 0]; 
epp=strtra2'' (-1) *eipp; 

%Reflection coefficients 
gamp=0 ; 

if (2*rs(m)+ctpp)~=0 
gamp=l/(2*rs(m)+ctpp); 
end 
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gamn=l/ (2*rs(m)*cos(thpp)+1); 

%Surface Currents 

Jxpp(i,m)=(epp(2)*gamp-epp(3)*gamn)*ctpp"2; 

Jypp(i,m)=(epp(2)*gamp+epp(3)*gamn)*ctpp; 

% Area integral for general case 

DD=Dq-Dp; 

expDo=exp(j * Do); 

expDp=exp(j*Dp); 

expDq=exp(j * Dq); 

% Special Case 1 
if abs(Dp) < Lt & abs(Dq)>= Lt 
sic=0.; 
for n=0:Nt 

sic=sic+(j*Dp)^n/fact(n)*(-Co/(n+1)+expDq*(Co*gf (n,-Dq) ) ) ; 
end 

Ic(i,m)=sic*2*Area(m)*expDo/j/Dq; 

% Special Case 2 

elseif abs(Dp) < Lt & abs(Dq) < Lt 
sic=0.; 
for n=0:Nt 

for nn=0:Nt 

sic=sic+ (j*Dp) ''n* (j*Dq) "'nn/fact (nn+n+2) *Co; 

end 

end 

Ic (i,in) =sic*2*Area (in) *expDo; 

% Special Case 3 

elseif abs(Dp) >= Lt & abs(Dq) < Lt 
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for n=0:Nt 


sic=sic+ (j*Dq) ''n/fact (n) *Co*gf (n+l, -Dp) / (n+1) ; 
end 

Ic(i,m)=sic*2*Area(m)*expDo*expDp; 

% Special Case 4 

elseif abs(Dp) >= Lt & abs(Dq) >= Lt & abs(DD) < Lt 
sic=0.; 
for n=0:Nt 

sic = sic+(j*DD)''n/fact (n) * (-Co*gf (n, Dq)+expDq*Co/(n+l) ) ; 
end 

Ic (i,in) =sic*2*Area (m) *expDo/j/Dq; 
else 

Ic (i,m) =2*Area (m) *expDo* (expDp*Co/Dp/DD-expDq*Co/Dq/DD-Co/Dp 
/Dq) ; 

end % end of special cases 

test 

%Scattered field in local rectangular coordinates 
Esxpp=Jxpp (i,in) *Ic (i,m) ; 

Esypp=Jypp(i,m)*Ic(i,m) ; 

Eszpp=0; 

%retranslation to global rectangular coordinates 
Espp(m, :) = [Esxpp Esypp Eszpp]; 

Esl=tra2'' (-1) *Espp (m, :) ' ; 

Es=tral"' (-1) *Esl; 

%retranslation to global spherical coordinates 
%seipp (m, :) =strtra'' (-1) *Es; 
seipp(m,1)=0; 
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seipp(m,2)=ct*cp*Es(1)+ct*sp*Es(2)-st*Es(3) ; 
seipp(m,3)=-sp*Es(1)+cp*Es(2); 

%suin over all triangles to get scattered field 
sumt=sunit+seipp (m, 2) ; 
suinp=sump+seipp (m, 3) ; 

end %end of illumination test 
end %end of triangle loop 
Eth=sumt; 

Sth(i)=10*logl0 (4*pi*abs (Eth) ''2+l.e-lO) ; 

Eph=sump; 

Sph(i)=10*logl0(4*pi*abs{Eph)^2+l.e-lO); 
end 

Sth{:)=max(Sth(:),-50) ; 

Sph(:)=max(Sph(:),-50); 

%plot res 
figure(2) 

plot(theta,Sth,theta,Sph,'—'),grid 

title{'solid: co-polarization dashed: cross-polarization') 
xlabel('monostatic angle, theta (deg)') 
ylabel('sigma/wavl^2 (dB)') 
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% P0TC02.M 

% PO scattering from triangles 
clear all 
format long 

% illumination flag: =0 external face only 
iflag=0; 

Lt=.001; 

Nt=2 ; 
wave=l; 
bk=2*pi/wave; 
rad=pi/180; 

% Incindent wave polarization 

Et=l+j*0; %TM-z 

Ep=0+j*0; %TE-z 

% Wave ampltidue at all vertices 

Co=l.; 

build 

ntria=nfaces; 
nvert=nverts; 


z=zpts; 

title('triangular model of the scattering surface') 

xlabel{'x') 

ylabel('y') 

zlabel('z') 

for i=linvert 

text(x(i)-max(x)/20,y(i)-max(x)/20,z(i),num2str(i)) 

end 

hold off 
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disp('hit return to continue calculation') 
pause 

% Define position vectors to vertices 
for i=l:nvert; 

r(i, :) = [X(i) y(i) z(i) ] ; 

end 

%Use res.m program to calculate resistivity 
%res 

% Get resistivity, edge vectors and normals from edge cross 

products 

for i=l:ntria 

rs(i)=0; 

%resistivity for a plate 

%rs(i)=rss*sqrt((xm(i)-.5)^2+(ym(i)-.5)^2)/377; 

% resistivity for a disk 

%rs (i) =rss*sqrt (xm(i) ''2+ym(i) ''2) / 377; 

% resistivity for a cap 
%rs (i)=rss*arc [±) / 311 ; 
ac(i)=l/(2*rs(i)+l) ; 

A=r(vind(i,2),:)-r(vind(i, 1) , ;) ; 

B=r(vind(i,3),;)-r(vind(i, 1) 

C=r(vind(i,3),:)-r(vind(i, 2) , ; 

N (i,:)=cross(A,B); 

% check for proper direction of normal vector 
if N(i,3)<0 

N(i, ;)=-N(i,:); 

end 

% Edge lengths for triangle "i" 
d(i,1)=norm(A); 
d(i, 2) =norm(B); 
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d(i, 3)=norm(C); 
ss=.5*sum(d(i,:)); 

Area(i)=sqrt(ss*(ss-d(i,1))*(ss-d(i, 2))*(ss-d(i,3))); 
Nn=norin (N (i, ;) ) ; 

N(i,:)=N(i,:)/Nn; 

%transform triangles into local coordinates 

xl=[l 0 0]; 

yl=[0 1 0 ] ; 

zl=[0 0 1] ; 

end 

% Pattern Loop 
it=0; 

load du.m; 
k=0; 

for ul=-l:du:1-du 

for vl=-l:du:1-du; 

it=it+l; 

Eth(it)=0; 

Eph(it)=0; 

%test to see if center of box is within visible region 
test={ (ul+du/2) "2+ (vl+du/2) ''2) ; 
if test<=l 
k=k+l; 

phi=atan2((vl+du/2),(ul+du/2)); 
cp=cos(phi); sp=sin(phi); 
theta=asin(sqrt(test)); 
st=sin(theta); ct=cos(theta); 

% Center of box is used so u and v are redefined using center 
of % box 

u=st*cp; v=st*sp; w=ct; 
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up=( 


ct*cp; vp=ct*sp; wp=-st; 

% Spherical coordinate system unit vectors 
R(it,:)=[u V w]; 

% Change Efield into global rectangular coordinates 

strtra(l,:)=[st*cp ct*cp -sp]; 

strtra(2,:)=[st*sp ct*sp cp]; 

strtraO, :) = [ct -st 0] ; 

sei=[0 Et Ep]; 

ei=strtra*sei'; 

% Begin loop over triangles 
for m=l:ntria 
t(m,;)=cross(zl, N(m, :) ) ; 
beta(m)=acos(N(m,3)); 

alpha(m)=atan2(t(m,2),t(m,1)+le-6)-pi/2; 
lc(k,m)=0; 

% Test to see if front face is illuminated 
ndotk=N(m,:)*R(it, :) '; 
if ndotk >= 0 & iflag==0 | iflag~=0 
%translation into local coordinates 
%first rotate about z axis 

tral(1,:)=[cos(alpha(m)) sin(alpha(m)) 0]; 

tral(2,;)=[-sin(alpha(m)) cos(alpha(m)) 0]; 

tral(3,:) = [0 0 1] ; 

% Rotation about yp axis 

tra2(1,:)=[cos(beta(m)) 0 -sin(beta(m))]; 

tra2(2,:)=[0 1 0]; 

tra2(3,:)=[ sin(beta(m)) 0 cos(beta(m))]; 
Dl=tral*R(it> 

D2=tra2*Dl; 
upp=D2(1); 
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vpp=D2(2); 
wpp=D2(3); 

% Find spherical angles in local coordinates 
stpp=sqrt (upp''2+vpp''2) *sign (wpp) ; 
ctpp=sqrt(l-stpp"2); 
thpp=acos (ctpp); 
phipp=atan2(vpp,upp+le-10); 

% Phase at the three vertices of triangle m; monostatic RCS 
needs % "2" 

Dp=2*bk* ( (x (vind (m, 1))-x(vind(m,3)))*u+.. . 

(y(vind(m,1))-y(vind(m,3)))*v+... 

(z (vind (m, 1) ) -z (vind (m, 3) ) ) *w) ; 

Dq=2*bk*((X(vind(m,2))-x(vind(m, 3)))*u+... 

(y(vind(m,2))-y(vind(m,3)))*v+. . . 

(z(vind(m,2))-z(vind(m, 3)))*w) ; 

Do=2*bk*(x(vind(m,3))*u+y(vind(m, 3) )*v+z(vind (m,3))*w) ; 

%incident field in local rectangular coordinates 

eip=tral*ei; 

eipp=tra2*eip; 

% translation of E-field into local spherical coordinates 
strtra2=[stpp*cos (phipp) ctpp*cos(phipp) -sin(phipp); 
stpp*sin(phipp) ctpp*sin(phipp) cos(phipp); ctpp -stpp 0]; 
epp=strtra2'' (-1) *eipp; 

% Area integral for general case 
DD=Dq-Dp; 
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expDo=exp (j*Do); 
expDp=exp(j*Dp); 
expDq=exp(j*Dq); 

% Special Case 1 
if abs(Dp) < Lt & abs(Dq)>= Lt 
sic=0.; 
for n=0:Nt 

sic=sic+ (j*Dp)^n/fact(n) *(-Co/(n+1) +expDq*(Co*gf(n,-Dq))) ; 
end 

Ic(k,m)=sic*2*Area(m)*expDo/j/Dq; 

% Special Case 2 

elseif abs(Dp) < Lt & abs(Dq) < Lt 
sic=0. ; 
for n=0;Nt 

for nn=0:Nt 

sic=sic+ (■j*Dp) ''n* (j*Dq) ''nn/fact (nn+n+2) *Co; 

end 

end 

Ic{k,m)=sic*2*Area(m)*expDo; 

% Special Case 3 

elseif abs(Dp) >= Lt & abs(Dq) < Lt 
sic=0.; 
for n=0:Nt 

sic=sic+ (j*Dq) '"n/fact (n) *Co*gf (n+.l,-Dp) / (n+1) ; 
end 

Ic (k,in) =sic*2*Area (m) *expDo*expDp; 

% Special Case 4 


83 






elseif abs(Dp) >= Lt & abs(Dq) >= Lt & abs(DD) < Lt 
sic=0.; 
for n=0:Nt 

sic = 

sic+ (j *DD) ''n/fact (n) * (-Co*gf (n, Dq) +expDq*Co/ (n+1) ) ; 
end 

Ic(k,m)=sic*2*Area(m)*expDo/j/Dq; 
else 

Ic(k,m)=2*Area(m)*expDo*(expDp*Co/Dp/DD-expDq*Co/Dq/DD-Co/Dp 
/Dq) ; 

end 

% end of special cases test 
end %end of illumination test 
end %end of triangle loop 

%sum over all triangles to get scattered field 
Eth(it)=Ic(k,:)*ac'*Et; 

Eph(it)=Ic(k,:)*ac'*Ep; 

% correction for rounding errors 
%Eth2(it)=round(Eth(it)*le4); 

%Eth(it)=Eth2(it)*le-4; 

end % end of test for visible region 

Sth(it)=10*logl0 {4*pi*abs {Eth(it) / (Et+le-10) ) ''2+l.e-lO) ; 

Sph (it) =10*logl0 (4*pi*abs (Eph (it) / (Ep+1 .e-10) ) ''2 + 1 .e-10) ; 
end % end of v loop 
end % end of u loop 
Sth(:)=max(Sth(:),-50); 

Sph(:)=max(Sph(:),-50); 

Sth=Sth'; 

Sph=Sph'; 

%plot res 
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u=-l;du;1-du; 
v=-l:du:1-du; 

T=ones(1,2/du); 
nrow=2/du; 
ncol=2/du; 
for k=l:ncol; 

k;l=nrow* (k-l)+l; 

k2=kl+nrow-l; 

X2(:,k)=(u(k)+du/2)*T'; 

Y2(:,k)=(du/2.+v) 

Z2(:,k)=Sth{kl:k2); 

end 

figure(2) 
orient tall 
surf(X2,Y2, Z2) 
view(-25, 55) 
colormap(hot) 
xlabel('u') 
ylabel('v') 
zlabel('RCS in Db’) 

Eth=Eth'; 

Eph=Eph'; 
rs=rs'; 

Ethr=real(Eth) ; 

Ethi=imag(Eth) ; 

Ephr=real(Eph) ; 

Ephi=imag(Eph) ; 

Ici=imag(Ic); 

save Ethr.m Ethr -ascii -double 
save Ethi.m Ethi -ascii -double 
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save Ephr.m Ephr -ascii 
save Ephi.m Ephi -ascii 
save rs.m rs -ascii 
save Icr.iti Ic -ascii 
save Ici.m Ici -ascii 
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% BUILD.M 

% Program to read transposed "buildnS" data and plot it in 
MATLAB 

% node (x,y,z) coordinates 
% legend: 
load xpts.m 
load ypts.m 
load zpts.m 
nverts=length(xpts); 

% node connection list 
load endl.m 
load end2.m 
nedges=length(endl) ; 

% face connection list 
load nodel.m 
load node2.m 
load nodes.m 
nfaces=length(nodes) ; 

% load vind array which gives the three nodes for each 
triangle 
for i=l:nfaces 

pts=[nodel(i) node2{i) nodeS(i)]; 
vind(i,:)=pts; 
end 

figured) 

% This section plots a mesh 
if icur==0 


for i=l:nfaces 

X= [xpts (vind (i, 1) ) xpts (’ 


d(i,2)) xpts(vind(i, 3) ) 




Y=[ypts(vind(i,1)) ypts(vind(i,2)) ypts(vind(i,3)) 
ypts(vind (i,1) ) ] ; 

Z=[zpts(vind(i,1)) zpts(vind(i,2)) zpts(vind(i,3)) 
zpts(vind(i,1))] ; 

plot3(X,Y,Z,'m') 
if i == 1 
axis equal 
view(60,45) 
hold on 
end 
end 
end 

hold off 
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% FACT.M 

% Subroutine for calculating the factorial of a number 
function f=fact(m) 



for n=l:m 

f=f*n; 

end 

end 






% GF.M 

function gf=gf(n,w) 

gf=(exp(j*w)-1)/(j*w) ; 
if n>=0 
for m=l:n 

gf=(exp(j*w)-n*gf)/(j*w) ; 
end 

end 
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% RES.M 

% Program to Calculate Resistivity on Each Triangle using a 
% Linear Taper 
for n=l:nfaces 

xm (n) = (r (max (vind (n, :) ) , 1) +r (min (vind (n, :) ) , 1) ) /2 ; 
ym (n) = (r (max (vind (n, :) ) , 2) +r (min (vind (n, : ) ) , 2) ) / 2; 
zm(n) = (r (max(vind(n,:)),3)+r(min(vind (n, :)),3))/2; 

% vector to midpoint of each triangle 
mn (n,:) = [xm(n) ym(n) zm(n)]; 

% angle from z axis to midpoint vector 
if zm(n)~=0 

omega (n) =pi/2-acos (zm (n) / (sqrt (xm(n) ''2+ym(n) ^2) + le-8) ) ; 
%arclength from z axis to midpoint of each triangle 
arc (n) =xm(n) ^2+ym (n) "'2; 
end 

%arclength from z axis to midpoint if plate is not flat 

arc(n)=omega(n)/(2*pi); 

end % end of triangle loop 

rss=377/max(abs(arc)); 

end 
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% SYNTH.M 

% Resistivity Synthesis-Monostatic 
clear all 
format long 

% illumination flag: =0 external face only 
iflag=0; 

Lt=.001; 

Nt=2 ; 
wave=l; 
bk=2*pi/wave; 
rad=pi/180; 

% Incindent wave polarization 
Et=l+j*0; %TM-z 

Ep=0+j*0; %TE-z 

% Load results from potcos program 
load Ethr.m; 
load Ethi.m; 

Ethii=Ethi*i; 

Eth=Ethr+Ethii; 
load Ephr.m; 
load Ephi.m; 

Ephii=Ephi*i; 

Eph=Ephr+Ephii; 
load rs.m; 

% Wave ampltidue at all vertices 

Co=l.; 

build 

ntria=nfaces; 
nvert=nverts; 
x=xpts; 
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y=ypts; 

z=zpts; 

title ('triangular model of the scattering surface') 

xlabel( 'X') 

ylabel('y') 

zlabel ('z') 

for i=linvert 

text (x(i)-max(x)/20,y(i)-max(x)/2 0,z(i),num2str(i)) 

end 

hold off 

dispt'hit return to continue calculation') 
pause 
% label 

% Define position vectors to vertices 
for i=linvert,• 

r(i,:)=[x(i) y(i) z(i)]; 

end 

% Get edge vectors and normals from edge cross products 
for i=lintria 

A=r(vind(i,2),i)-r(vind(i,1),i); 

B=r(vind(i,3),i)-r(vind{i,1),i); 

C=r(vind(i, 3) , i)-r(vind(i,2),i); 

N(i,I)=cross(A,B); 

% checlc for proper direction of normal vector 
if N(i,3)<0 

N(i, I )=-N(i, I ) ; 

end 

% Edge lengths for triangle "i" 
d(i,l)=norm(A); 
d(i, 2)=norm(B); 
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d(i, 3)=norm(C); 
ss=.5*sum(d(i,;)); 

Area(i)=sqrt(ss*(ss-d(i,1))*(ss-d(i,2))*(ss-d (i,3))); 
Nn=norni (N (i, : ) ) ; 

N(i, :)=N (i,:)/Nn; 

%transform triangles into local coordinates 

xl=[l 00]; 

yl=[0 1 0]; 

zl=[0 0 1]; 

end 

% Pattern Loop 
it=0; 
k= On¬ 
load du.m; 
for ul=-l:du:1-du 

for vl=-l:du:1-du; 

it=it+l; 

%test to see if center of box is within visible region 
test= ( (ul+du/2) ^^2+ (vl+du/2) ''2) ; 
if test<=l 
k=k+l; 

phi=atan2((vl+du/2),(ul+du/2)); 
cp=cos(phi); sp=sin(phi); 
theta=asin(sqrt(test)); 
st=sin(theta); ct=cos(theta); 

% Center of box is used so u and v are redefined using center 
of % box 

u(k)=st*cp; v(k)=st*sp; w=ct; 
up=ct*cp; vp=ct*sp; wp=-st; 

% Spherical coordinate system unit vectors 
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R{it,:)=[u(k) v(k) w]; 

% Change Efield into global rectangular coordinates 

strtra(1,;)=[st^cp ct*cp -sp]; 

strtra(2,:)=[st*sp ct*sp cp]; 

strtraO, :) = [ct -st 0] ; 

sei=[0 Et Ep]; 

ei=strtra*sei'; 

% Begin loop over triangles 
suint=0; 
sump=0; 
for m=l;ntria 
t (m, : ) =cross (zl, N (m, : ) ) ; 
beta(m)=acos(N(m,3)); 

alpha(m)=atan2(t(m,2),t(m,l)+le-6)-pi/2; 

% Test to see if front face is illuminated 
ndotk=N(m,:)*R(it,:)'; 
if ndotk >= 0 & iflag==0 | iflag~=0 
%translation into local coordinates 
%first rotate about z axis 

tral (1,:) = [cos(alpha(m)) sin(alpha(m)) 0]; 

tral(2,:) = [-sin(alpha(m)) cos(alpha(m)) 0] ; 

tral (3, :) = [0 0 1]; 

% Rotation about yp axis 

tra2(1,:)=[cos(beta(m)) 0 -sin(beta(m))]; 
tra2 (2,:) = [0 1 0]; 

tra2(3,:)=[ sin(beta(m)) 0 cos(beta(m))]; 
Dl=tral*R(it,:)'; 

D2=tra2*Dl; 
upp=D2(1); 
vpp=D2(2); 
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wpp=D2(3); 

% Find spherical angles in local coordinates 
stpp=sqrt (upp"'2+vpp''2) *sign (wpp) ; 
ctpp=sqrt(l-stpp^2); 
thpp=acos(ctpp); 
phipp=atan2(vpp,upp+le-10); 

% Phase at the three vertices of triangle m; monostatic RCS 
needs % "2" 

Dp=2*bk* ( (x (vind (m, 1) ) -x (vind (m, 3) ) ) *u (k) + . . . 

(y(vind(m,1))-y(vind(m,3)))*v{k)+... 

(z(vind{m,1))-z(vind(m,3)))*w); 

Dq=2*bk* ( (x (vind (m, 2) ) -x (vind (m, 3) ) ) *u (k) + . . . 

(y(vind(m,2))-y(vind(m,3)))*v(k)+... 

(z (vind (m, 2) ) -z (vind (m, 3) ) ) *w) ; 

Do=2*bk* (X (vind (m, 3) ) *u (k) +y (vind (m, 3) ) *v (k) +z (vind (m, 3) ) *w) ; 

%incident field in local rectangular coordinates 

eip=tral*ei; 

eipp=tra2*eip; 

% translation of E-field into local spherical coordinates 
strtra2=[stpp*cos(phipp) ctpp*cos(phipp) -sin(phipp) ; 
stpp*sin(phipp) ctpp*sin(phipp) cos(phipp); ctpp -stpp 0); 
epp=strtra2'' (-1) *eipp; 

% Area integral for general case 
DD=Dq-Dp; 
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expDo=exp (j *Do); 
expDp=exp(j*Dp); 
expDq=exp (j *Dq); 

% Special Case 1 
if abs(Dp) < Lt & abs(Dq)>= Lt 
sic=0.; 
for n=0;Nt 

sic=sic+ {j*Dp) ''n/fact (n) * (-Co/ (n+1) +expDq* (Co*gf (n, -Dq) ) ) ; 
end 

Ic(k,m)=sic*2*Area(m)*expDo/j/Dq; 

% Special Case 2 

elseif abs(Dp) < Lt & abs(Dq) < Lt 
sic=0.; 
for n=0;Nt 

for nn=0:Nt 

sic=sic+ (j *Dp)^n*( j *Dq)^nn/fact(nn+n+2)*Co; 

end 

end 

Ic (k,iti) =sic*2*Area (m) *expDo; 

% Special Case 3 

elseif abs(Dp) >= Lt & abs(Dq) < Lt 
sic=0.; 
for n=0:Nt 

sic=sic+ (j*Dq) ''n/fact (n) *Co*gf (n+1, -Dp) / (n+1) ; 
end 

Ic(k,m)=sic*2*Area(m)*expDo*expDp; 

% Special Case 4 
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elseif abs(Dp) >= Lt & abs(Dq) >= Lt & abs(DD) < Lt 
sic=0.; 
for n=0:Nt 

sic = sic+(j*DD)''n/fact (n) * (-Co*gf (n, Dq)+expDq*Co/(n+1) ) ; 
end 

Ic(k,m)=sic*2*Area(m)*expDo/j/Dq; 
else 

Ic (k,m) =2*Area (iti) *expDo* (expDp*Co/Dp/DD-expDq*Co/Dq/DD-Co/Dp 
/Dq) ; 

end 

% end of special cases test 
end %end of illumination test 
end %end of triangle loop 
Eth2(k)=Eth(it); 

end % end of test for visible region 

end % end of v loop 

end % end of u loop 

ac=pinv(Ic)*Eth2'; 

for k=l:ntria 

if imag(ac(k))<=.0001 
ac(k)=real(ac(k)); 
end 

sr(k)=(l/ac(k)-l)*.5; 

% correct round-off 
s (k)=round(sr(k)*le4) ; 
sr{k)=s(k)*le-4; 
end 

k=l:ntria; 
figure(2) 

plot(k,rs,k,sr,'-'); 
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title{'Comparison of Original and Synthesized Resistivity') 
xlabel('Triangle Number') 
ylabel('Normalized Resistivity') 
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% POTESS.M 

% PO scattering from triangles 
clear all 

% illumination flag; =0 external face only 
iflag=l; 

Lt=.001; 

Nt=2 ; 
wave=l; 
bk=2*pi/wave; 
rad=pi/180; 

% Incindent wave polarization 
Et=l+j*0; %TM-z 

Ep=0+j*0; %TE-z 


phr=phi*rad; 

cp=cos(phr);sp=sin(phr); 

% Wave ampltidue at all vertices 

Co=l.; 

build 

ntria=nfaces; 
nvert=nverts; 
x=xpts; 
y=ypts; 
z=zpts; 

title('triangular model of the scattering surface') 

xlabel('X') 

ylabel('y') 

zlabel('z') 

for i=l;nvert 

text(x(i)-max(x)/20,y(i)-max(x)/20,z(i), num2str ( i)) 
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end 

hold off 

dispC'hit return to continue calculation') 
pause 

% Define position vectors to vertices 
for i=l;nvert; 

r(i,:)=[x(i) y{i) z(i)]; 

end 

%Use synthesized resistivity 

%load srO.m % 0 resistivity 

%load srl.m % 377 ohm/sq resistivity 

load srl.m % linear resistivity 

rs=srl; 

% edge vectors and normals from edge cross products 
for i=l:ntria 

A=r(vind(i,2),:)-r(vind(i, 1), ;) ; 

B=r(vind(i,3),:)-r(vind(i,1),:); 

C=r(vind(i,3),:)-r(vind(i,2), :) ; 

N(i,;)=cross(A,B); 

% check for proper direction of normal vector 
if N(i,3)<0 

N(i,;)=-N(i,:); 

end 


% Edge lengths for triangle "i" 
d(l,1)=norm(A) ; 
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d(i,2)=norm(B); 
d(i,3)=norm(C); 
ss=.5*sum(d(i, :) ) ; 

Area(i)=sqrt(ss*(ss-d(i,1))*(ss-d(i, 2) )*(ss-d(i,3))) 
Nn=norm(N(i,:)); 

N(i,:)=N(i,:)/Nn; 

%transforin triangles into local coordinates 

xl=[l 0 0]; 

yl=[0 1 0]; 

zl=[0 0 1]; 

end 

% Pattern Loop 
start=0; 
stop=180; 
del=l; 

it=floor((stop-start)/del) +1; 
for i=l:it 

theta (i)=start+(i-1)*del ; 
thr=theta(i)*rad; 
st=sin(thr); ct=cos(thr); 
u=st*cp; v=st*sp; w=ct; 
up=ct*cp; vp=ct*sp; wp=-st; 

% Spherical coordinate system unit vectors 
R(i, :)=[u V w]; 

% Change Efield into global rectangular coordinates 

strtra(1,:)=[st*cp ct*cp -sp]; 

strtra(2, : ) = [st*sp ct*sp cp]; 

strtra(3,:)=[ct -st 0]; 

sei= [0 Et Ep]; 
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ei=strtra*sei'; 

% correct orientation for normal after 90deg 
if theta(i)==91 
N=-N; 

end 

% Begin loop over triangles 
sumt=0; 
sui[ip=0; 
for m=l:ntria 
t (m, : ) =cross (zl, N (m, :) ) ; 
beta(m)=acos(N(m,3)); 

alpha(m)=atan2(t(m,2),t(ra,1)+le-6)-pi/2; 

% Test to see if front face is illuminated 
ndotk=N(m,;)*R(i,:)'; 
if ndotk >= 0 & iflag==0 | iflag~=0 
%translation into local coordinates 
%first rotate about z axis 

tral(1,:) = [cos (alpha(m)) sin(alpha(m)) 0]; 
tral (2, :) = [-sin(alpha(m)) cos(alpha(m)) 0] ; 
tral(3,:)=[0 01]; 

% Rotation about yp axis 

tra2(1,:)=[cos(beta(m)) 0 -sin(beta(m))]; 
tra2(2,:)=[0 1 0]; 

tra2(3,:)=[ sin(beta(m)) 0 cos(beta(m))]; 
Dl=tral*R(i,:)'; 

D2=tra2*Dl; 
upp=D2(1); 
vpp=D2(2); 
wpp=D2(3); 

% Find spherical angles in local coordinates 
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stpp=sqrt (upp"'2+vpp''2) *sign (wpp) ; 
ctpp=sqrt (l-stpp''2) ; 
thpp=acos(ctpp); 
phipp=atan2(vpp,upp+le-10); 

% Phase at the three vertices of triangle in; monostatic RCS 
needs % "2" 

Dp=2*bk*((x(vind(m,1))-x(vind (m, 3)))*u+ . . . 

(y(vind(m,1))-y(vind(m,3)))*v+... 

(z(vind(m,1))-z(vind(m,3)))*w); 

Dq=2*bk*((X(vind(m,2))-x(vind(m,3)))*u+ . . . 

(y (vind (m, 2) ) -y (vind (m, 3) ) ) *v+ . . . 

(z (vind (m, 2) ) -z (vind (m, 3) ) ) *w) ; 

Do=2*bk* (x(vind(m,3))*u+y(vind(m, 3))*v+z(vind (m, 3))*w); 
%incident field in local rectangular coordinates 
eip=tral*ei; 
eipp=tra2*eip; 

% translation of E-field into local spherical coordinates 
strtra2=[stpp*cos(phipp) ctpp*cos(phipp) -sin(phipp); 
stpp*sin(phipp) ctpp*sin(phipp) cos(phipp); ctpp -stpp 0]; 
epp=strtra2'' (-1) *eipp; 

%Reflection coefficients 
gainp=0; 

if (2*rs(m)+ctpp)~=0 
gamp=l/(2*rs(m)+ctpp); 
end 

gainn=l/ (2*rs (m) *cos (thpp) +1) ; 
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^Surface Currents 

Jxpp (i,m) = (epp (2) *gamp-epp (3) *gaiim) *ctpp''2; 

Jypp (i, m) = (epp (2) *gamp+epp (3) *gainn) *ctpp; 

% Area integral for general case 

DD=Dq-Dp; 

expDo=exp (j*Do) ; 

expDp=exp(j*Dp); 

expDq=exp(j*Dq); 

% Special Case 1 

if abs(Dp) < Lt & abs(Dq)>= Lt 

%disp('case 1') 

lease(i,m)=l; 
sic=0.; 
for n=0:Nt 

sic=sic+ (j *Dp)^n/fact(n)*(-Co/(n+1)+expDq*(Co*gf(n, -Dq) ) ); 
end 

Ic (i,m)=sic*2*Area(m)*expDo/j/Dq; 

% Special Case 2 

elseif abs(Dp) < Lt & abs(Dq) < Lt 
%disp('case 2') 

lease (i,m)=2; 
sic=0.; 
for n=0:Nt 

for nn=0;Nt 

sic=slc+(j*Dp)^n*(j*Dq)^nn/fact(nn+n+2)*Co; 
end 
end 

Ic(i,m)=sic*2*Area(m)*expDo; 
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% Special Case 3 

elseif abs(Dp) >= Lt & abs(Dq) < Lt 
%disp('case 3') 

lease(i,m)=3; 
sic=0.; 
for n=0:Nt 

sic=sic+(j*Dq)^n/fact(n)*Co*gf(n+1,-Dp)/(n+1); 
end 

Ic(i,m)=sic*2*Area(m)*expDo*expDp; 

% Special Case 4 

elseif abs(Dp) >= Lt & abs(Dq) >= Lt & abs(DD) < Lt 
%disp('case 4') 

lease(i,m)=4; 
sic=0.; 
for n=0:Nt 

sic = sic+(j*DD)^n/fact(n)*(-Co*gf(n,Dq)+expDq*Co/(n+1) ) ; 
end 

Ic (i,in) =sic*2*Area (m) *expDo/j/Dq; 
else 

lease(i,m)=0; 

Ic(i,m)=2*Area(m)*expDo*(expDp*Co/Dp/DD-expDq*Co/Dq/DD-Co/Dp - 
/Dq) ; 

end 

% end of special cases test 

%Scattered field in local rectangular coordinates 
Esxpp=Jxpp (i,in) *Ic(i,m) ; 

Esypp=Jypp(i,m)*Ic(i,m); 

Eszpp=0; 

%retranslation to global rectangular coordinates 
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Espp(m,:)=[Esxpp Esypp Eszpp]; 

Esl=tra2^(-1)*Espp(m,:)'; 

Es=tral'' (-1) *Esl; 

%retranslation to global spherical coordinates 
%seipp (m, : ) =strtra'' (-1) *Es; 
seipp(m,1)=0; 

seipp(m,2)=ct*cp*Es(1)+ct*sp*Es(2)-st*Es(3); 
seipp(m,3)=-sp*Es(1)+cp*Es(2); 

%sum over all triangles to get scattered field 

sumt=sumt+seipp(m,2); 

smnp=suinp+seipp (m, 3) ; 

end %end of illumination test 

end %end of triangle loop 

Eth=sumt; 

Sth ( i) =10*logic (4*pi*abs (Eth) ''2+1. e-10) ; 
Eph=sump; 

Sph (i) =10*logl0 (4*pi*abs (Eph) ''2+1.e-10) ; 
end 

Sth(:)=max(Sth(:),-50); 

Sph(:)=max(Sph(:),-50); 

%load rcsO.m 
%load rcsl.m 
load rcsl.m 
%plot res 
figure (2) 

plot(theta,rcsl,theta,Sth,'—'),grid 
title('solid:original dashed:synthesized') 
xlabel('monostatic angle, theta (deg)') 
ylabel ( ' sigma/wavl''2 (dB) ' ) 
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